arXiv:1507.05590v2 [physics.atom-ph] 25 Sep 2015 


First-principles nonequilibrium Green’s function approach to transient photoabsorption: 

Application to atoms 


E. Perfetto/’^ A.-M. Uimonen,^ R. van Leeuwen,^ and G. Stefanucci^’^ 

^Dipartimento di Fisica, Universitd di Roma Tor Vergata, Via della Ricerca Scientifica 1, 00133 Rome, Italy 
^INFN, Laboratori Nazionali di Frascati, Via E. Fermi 40, 00044 Frascati, Italy 
^Department of Physics, Nanoscience Center, FIN 40014, University of Jyvdskyld, Jyvdskyld, Finland 
Department of Physics, Nanoscience Center, FIN 40014, University of Jyvdskyld, 

Jyvdskyld, Finland; European Theoretical Spectroscopy Facility (ETSF) 

^Dipartimento di Fisica, Universitd di Roma Tor Vergata, Via della Ricerca Scientifica 1, 

00133 Rome, Italy; European Theoretical Spectroscopy Facility (ETSF) 

We put forward a first-principle NonEquilibrium Green’s Function (NEGF) approach to calculate the transient 
photoabsorption spectrum of optically thin systems. The method can deal with pump fields of arbitrary strength, 
frequency and duration as well as for overlapping and nonoverlapping pump and probe pulses. The electron- 
electron repulsion is accounted for by the correlation self-energy, and the resulting numerical scheme deals 
with matrices that scale quadratically with the system size. Two recent experiments, the first on helium and 
the second on krypton, are addressed. For the first experiment we explain the bending of the Autler-Townes 
absorption peaks with increasing the pump-probe delay r, and relate the bending to the thickness and density 
of the gas. For the second experiment we find that sizable spectral structures of the pump-generated admixture 
of Kr ions are fingerprints of dynamical correlation effects, and hence they cannot be reproduced by time-local 
self-energy approximations. Remarkably, the NEGF approach also captures the retardation of the absorption 
onset of Kr^^ with respect to Kr^^ as a function of r. 
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I. INTRODUCTION 

Transient photoabsorption (TPA) spectroscopy has today 
become a popular technique to investigate the ultrafast dy¬ 
namics of electrons and nuclei in atoms, molecules and 
nanostructures.^"^ A reliable physical interpretation of the 
TPA spectrum is inescapably linked to a reliable calculation 
of the probe-induced polarization in the pump-driven sys¬ 
tem. State-of-the-art calculations are based on the Configu¬ 
ration Interaction (Cl) expansion of the time-evolved many- 
electron state. The time-dependent Cl coefficients are ei¬ 
ther varied with respect to the probe field^"^^ or used to 
construct the dipole response function from the Lehmann 
representation, the latter approach being applicable only 
provided that the dressed pump and probe fields do not 
overlap.However, the size of the arrays in Cl calculations 
scale exponentially with the number of basis functions and 
the time-step to achieve convergence is typically much smaller 
than the time-step used in statistical approaches. Due to these 
numerical limitations the Cl approach is confined to the study 
of rather small systems. 

One possible statistical approach to TPA spectroscopy is 
Time Dependent Density Eunctional Theory (TDDET).^^ In 
TDDET are the occupied single-particle wavefunctions that 
are propagated in time and since they scale linearly with the 
number of basis functions TDDET is suitable to study much 
larger systems than those accessible by CL In the frame¬ 
work of the Adiabatic Local Density Approximation (ALDA) 
TDDET has been recently and successfully applied to the 
study of TPA in small- and medium-sized moleculesas 
well as to monitor the vibronic-mediated charge transfer in 
donor-acceptor complexes.Still, ALDA functionals have 
drawbacks that could compromise the description of pho¬ 


toabsorption spectra even in equilibrium systems. Eor ex¬ 
ample, ALDA misses correlation-induced spectral features 
like double-excitations^^’^^ and long-range charge-transfer 
excitations it also provides a poor description of the 
energy-level alignment in metal/molecule interfaces'^and 
of the Coulomb blockade phenomenon.In this work we 
discover yet another correlation effect missed by ALDA. 

An alternative statistical approach to TDDET is the many- 
body diagrammatic theory. Here the building blocks of 
the formalism are the NonEquilibrium Green’s Eunctions 
(NEGE),^^"^^ and correlation effects are included by a proper 
selection of self-energy diagrams. Double-excitations and 
other properties missed by ALDA are within reach of dia¬ 
grammatic theory already with basic self-energies. Recently, 
it has been shown that the TPA spectrum follows from the 
solution of a nonequilibrium Bethe-Salpeter equation (BSE) 
provided that the time-scale of the pump-induced electron 
dynamics is much longer than the life-time of the dressed 
probe field.^^ In general, however, for pump fields of arbitrary 
strength, frequency and duration and/or for overlapping pump 
and probe pulses, the nonequilibrium BSE is inadequate and 
the full time-propagation of the NEGE is unavoidable. Never¬ 
theless, as the size of the arrays in NEGE calculations scales 
quadratically with the number of basis functions, this formal¬ 
ism too allows for extending the range of Cl accessible sys¬ 
tems. 

In this work we formulate a general first-principle NEGE 
scheme to TPA and apply it to reproduce the transient spectra 
of a thick helium gas^^ and of a krypton gas.^"^ Eor helium we 
address the exponential damping of the probe-induced dipole, 
and relate it to the thickness and density of the gas. We then 
provide the explanation of the bending of the Autler-Townes 
absorption peaks with increasing the pump-probe delay. We 
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also propose a useful formula for fitting the experimental TPA 
spectra. The krypton gas constitutes a more severe test for 
NEGF due to non-trivial correlation effects. We find that for 
a proper description of the (pump-generated) evolving admix¬ 
ture of Kr ions the self-energy should have memory. Static (or 
adiabatic) approximations like the Hartree-Fock or the Marko¬ 
vian approximations perform rather poorly as only the spec¬ 
trum of Kr^+ is visible. Instead, the TPA spectrum calculated 
using the (memory-dependent) second-Bom self-energy con¬ 
tains absorption lines attributable to excitations in the Kr^+ 
and Kr^+ ions. Remarkably, we are also able to reproduce the 
femtosecond retardation of the absorption onset of Kr^+ with 
respect to Kr^+ as a function of the pump-probe delay. 

The paper is organized as follows. In Section II we relate 
the TPA spectrum to the microscopic quantum mechanical av¬ 
erage of the transverse, probe-induced dipole-wave propagat¬ 
ing toward the detector. In Section III we put forward the 
NEGF approach to TPA and provide the explicit expression 
of the self-energy approximations implemented in this work. 
The helium gas is studied in Section IV. In Section V we 
extend the NEGF approach to deal with pump-induced ion¬ 
ization processes and then apply the extended approach to the 
study of the krypton gas in Section VI. Summary and conclu¬ 
sions are drawn in Section VII. 


II. TRANSIENT PHOTOABSORPTION SPECTRUM 

We consider a gas of atoms or molecules perturbed by a 
strong time-dependent transverse electric field E(rf) (pump) 
propagating along the unit vector N and a feeble time- 
dependent transverse electric field e(rt) (probe) propagating 
along the unit vector n, see Fig. 1. The direction of propaga¬ 
tion N / n and in experiments the photodetector is positioned 
along the probe beam-line. Fet Sp{Tt) be the component of 
the total electric field (external plus induced) propagating to¬ 
ward the detector and Sp{t) be the value of Sp{r t) at the de¬ 
tector surface. Then the transmitted energy measured by the 
detector is 

with S the surface of the sample (assumed to be smaller than 
the laser beam cross section). Here and in the following we 
use the convention that quantities with the tilde symbol on 
top denote the Fourier transform of the corresponding time- 
dependent quantities. Replacing Sp in Eq. (1) with the exter¬ 
nal probe field e we obtain the energy Ej of the incident probe 
beam. The absorbed energy Ea is defined according to 

Ea = Ej — Et 

By means of a spectrometer it is possible to measure the ab¬ 
sorbed energy per unit frequency, i.e., 

6(w) ='5^ (|eH|2 - |£p(w)|2^ . (3) 



FIG. 1: (Color online) Schematic illustration of a pump-probe exper¬ 
iment. 

The quantity &{uj) is the transient photoabsorption (TPA) 
spectrum that we are interested in. 

As both E and e are transverse fields we have E(rt) = 
E(X t) with A = N • r and e(r t) = e{x t) with x = n • r. By 
definition Sp too depends only on x, i.e., Sp{rt) = Sp{xt), 
since it is a transverse field propagating along n. Without loss 
of generality we choose the coordinates of the boundaries of 
the sample in x = 0 and x = i, i being the thickness of the 
gas, see Fig. 1. From Maxwell equations the relation between 
the total field and the incident probe field is 

Sp{xt) = e{xt)[ dx'dp{x't) (4) 
c at Jo 

where dp{xt) is the component of the probe-induced dipole 
density propagating toward the detector.Notice that dp{x t) 
vanishes for x > i; hence Sp{i t) and Sp{t) (the electric field 
at the detector surface) differ only by a time-shift. This time- 
shift is completely irrelevant for the calculation of the spec¬ 
trum, so we can either use Sp{iuj) or Sp{uj) in Eq. (3). 

Equation (4) connects the experimental outcome &{uj) 
to a quantum-mechanical average. Fet us define dp{xt) 
more rigorously. We denote by |T^(rt)) the many-body 
state of an atom located in a volume element around r at 
time t when both pump and probe fields are present; sim¬ 
ilarly |T^p(rf)) is the many-body state of the same atom 
when only the pump is present (probe-free state). The av¬ 
erage of the atomic dipole operator d(r) over these two 
states is d*^^*^(rt) = (T^(rt)|d(r)|T^(rt)) and dp^^(rt) = 
(T^p(rt)|d(r)|T^p(r t)) respectively. The atomic probe- 
induced dipole is therefore 

(r t) = (r t) - (r t) . (5) 

For isotropic systems the probe-free dp*^ (r t) = dp*^ {X t) 
is a transverse field propagating along N. Instead the atomic 
dipole is a transverse field propagating along all 

possible directions kn AN with k and K real numbers. 
To first order in |e| the component propagating along N 
(k = 0) is exactly dp^\xt). It is reasonable to expect 
that for /c / 0 and for small |e| the dominant component of 

d(at)(j.^) 

is the one propagating along the same direction n 
of the external probe. In this case the probe-induced atomic 
dipole dp^^^ (r t) = dp^^^ {x t) is a transverse field propagat¬ 
ing along n and the function dp appearing in Eq. (4) is simply 
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^(a t)d^at) 

, being the density of the gas. We can therefore 

calculate dp by performing a time-propagation with pump and 
probe, another time-propagation with only the pump and then 
subtracting the resulting dipoles. Below we derive the basic 
equations to perform these calculations. 

For an atom with unperturbed Hamiltonian H (r) the evolu¬ 
tion of the state | ^(r t)) is governed by the Schrodinger equa¬ 
tion 


i-|5'(rt)) = [H{r) + £{rt) • d(r)J |5'(rt)) (6) 


where S is the total electric field. We choose the pump and 
probe propagation directions such that N • n 1. Then 
X = r- N:^r-n = x inside the sample and every r- 
dependent quantity can be approximated with an x-only de¬ 
pendent quantity. In particular the total field can be approxi¬ 
mated as 


S{xt) = E(xt) + e{xt) 


c dt 


t). 


(7) 

For optically thin samples the dipole oscillations in two points 
xi and X 2 are time-shifted by an amount \xi — X 2 \/c which 
is much smaller than the inverse of the typical dipole fre¬ 
quencies. Therefore t) is weakly depen¬ 

dent on x' inside the gas and we can approximate the integral 
in Eq. (7) with xd^^^\t). The mathematical simplification 
brought about by this approximation is that every atom can be 
evolved separately. Consider the state |T^(t)) = \'^{iyzt)) 
of an atom at the interface inx = i and let ^ = H{iy z) and 
d = d{ty z)hQ the atomic Hamiltonian and dipole operator. 
Then for r = {iy z) Eq. (6) reads 


4 mm 


H + e{t)-d 


( 8 ) 


with 


The numerical solution of Eqs. (8-10) is impractical for 
heavy atoms or moderate size molecules. In the next section 
we describe a statistical approach which avoids solving the 
Schrodinger equation for the many-electron system. This ap¬ 
proach is based on the Nonequilibrium Green’s Functions^^"^^ 
(NEGF) combined with the Generalized Kadanoff-Baym 
Ansatz^^ (GKBA), and it has been successfully applied to 
the electron gas,^^’^^ two-band model semiconductors^^’^^""^^ 
and more recently bulk Silicon,Hubbard chains"^^""^^ and 
donor-acceptor junctions."^^ Here we extend it to perform first- 
principles simulations of laser-driven quantum systems. 


III. NEGF APPROACH 


We work in the formalism of second quantization and de¬ 
note by Qcr (c]^) the annihilation (creation) operator for an 
electron in the orbital ipi{r) with spin a =t,i- Without loss 
of generality the basis is taken orthonormal. The unper¬ 
turbed Hamiltonian reads 

H = 2 "^ijmn^l^C^j^/C^cr'Cna ( 13 ) 

ij ijmn 


with hij = /dr (/?*(r)[— ^ + 14(r) the one-electron 

integrals with nuclear potential 14 and 


Vijmn = I dr dr 






(14) 


the four-index Coulomb integrals. Similarly the dipole opera¬ 
tor reads 


d = (15) 

ij 


£{t) = 

c dt 

and 


with dij = fdr(p*(r)r(pj(r) the matrix elements of the 
dipole vector. The time-dependent average of the atomic 
dipole, see Eq. (10), is therefore 


= {^{t)\d\'^{t)). (10) 

Equations (8-10) form a close system of three coupled equa¬ 
tions. The same equations could be solved for e = 0 to obtain 
the probe-free atomic dipole dp^\t) and hence the probe- 
induced dipole dl^^\t) in accordance with Eq. (5). Having 
we can calculate the probe-induced dipole density 
dp{£t) = and subsequently the component of 

the electric field propagating toward the detector from Eq. (4), 
i.e., 

Ptt f d 

Sp{it) = e{it) + —^—dp{it). (11) 

Fourier transforming Sp{it) and e{it) the TPA spectrum 
&{uj) in Eq. (3) follows: 


ij 

(j 

Let us introduce the evolution operator Z7(t, to) a time 
to when the system is unperturbed to an arbitrary time t. Then 
|T^(t)) = t7(t, to)|T^o) with |T^o) the state of the unperturbed 
system. To distinguish the Heisenberg from the Schrodinger 
picture we insert a subscript “77” to an operator 0(t) in the 
Schrodinger picture, Onif) = 7/(to,t)0(t)7/(t,to). The 
lesser (G^) and greater (G^) Green’s functions are defined 
according to 

(17a) 

Gf^tX) = (17b) 


©(cc) = — 2Im 00 (£uj) ' dp(£uj) 


27 r£ 

Sc 


; dp{£oo) 


( 12 ) 


and have the property 


GfAt,t') = -[GtX,t )]4 


(18) 
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implying that from with times t > t' ov t < t' wt can 
reconstruct the entire G^. In this work we consider spin- 
compensated systems with no spin-orbit interaction; hence 
G^ takes the same value for a =t, The NEGF formalism, 
however, is not limited to this case and it can easily be general¬ 
ized to Hamiltonians that are not diagonal in spin space. From 
at equal times we can calculate the time-dependent aver¬ 
age of any one-body operator; in particular the atomic dipole 
in Eq. (16) reads = -2i J2ij t). 

The lesser and greater Green’s functions satisfy nonlinear 
integro-differential equations known as the Kadanoff-Baym 
equations At present the numerical solution of the 

KBE for inhomogeneous systems is possible only for moder¬ 
ate size basis sets,"^^"^^ and still the CPU time is of the order of 
a few days for a propagation of ^ 10^ time steps and ^10^ 
basis functions. Considering that a TPA spectrum typically 
requires 10^ time steps for every delay between the pump and 
probe pulses, the KBE are not feasible in this context. A way 
to drastically reduce the computational effort consists in mak¬ 
ing the GKBA^^ 

G^{t,t') = ^ p{t)G^{t,t'), (19a) 

G^{t,t') = G^{t,t')p{t') — p{t)G^{t,t'). (19b) 

In these equations G^(t,t') = [G^{t',t)]^ is the retarded 
Green’s function and p{t) = —iG^{t,t) = 1 — p{t) = 
1 — iG^ (t, t) is the one-particle density matrix (the quantity of 
interest for the calculation of the atomic dipole). The GKBA 
is exact in the Hartree-Fock (HE) approximation and it is ex¬ 
pected to be accurate in systems with well-defined quasipar¬ 
ticles, see Ref. 48 for a recent discussion. From the KBE we 
can easily derive, see Appendix A, the following equation of 
motion for p (in matrix form) 

+ [hHF{t),p{t)] = il{t) - H.C., (20) 

where the HE single-particle Hamiltonian reads 

~ hij • d^jf + ^ ^ "^imnjPnm (t ), (21) 

mn 

with 

^imnj = (22) 

and the collision integral reads (in matrix form) 

I{t) = [ + . 

J —CO 

(23) 


The correlation self-energy appearing in Eq. (23) is a 
functional of G^ and G^ which, in turn, are functionals of 
p through the GKBA. Thus, Eq. (20) is a closed nonlinear 
integro-differential equation for p once we specify the func¬ 
tional form of the self-energy and the retarded Green’s func¬ 
tion. 

In this work we present results obtained using the second- 
Born (2B) self-energy 



FIG. 2: Diagrammatic representation of the 2B self-energy. Wiggly 
lines denote the Coulomb interaction v. 


whose diagrammatic representation is shown in Fig. 2. For 
the retarded Green’s function we consider 


G^{t,t') 


-i0{t -1') T 


e 


i dthupi^t) 


(25) 


where T is the time-ordering operator. In Ref. 48 we dis¬ 
cussed how to include correlation effects (beyond HE) in G^ 
and showed their importance in quantum transport calcula¬ 
tions. For the finite systems analyzed here, however, we found 
that the HE G^ of Eq. (25) is accurate enough. Similar find¬ 
ings were recently found in strongly correlated systems as 
weiF56,57 addition to the simplicity of Eq. (25), the use of a 
HE G^ within the NEGF-fGKBA scheme guarantees the sat¬ 
isfaction of all basic conservation laws."^^ 

We emphasize that the equation of motion for p is an equa¬ 
tion with memory since the evaluation of the collision integral 
at time t involves the density matrix at all times t' < t. Us¬ 
ing the GKBA expression for the lesser and greater Green’s 
function the collision integral can be rewritten as 


Iik{t) = ^^VirpnWmqsj [ dt' {G^ {t, t') p{t'{G^ {t, t') p{t'{p{t')G^ {t' ,t)) {p{t')G^ {t' ,t)) 

n.m. A 'J —OO 


pq 

sr 


+ {G^{t,t')p{t'))^^ {G^{t,t')p{t')) {p{t')G^{t',t)){p{t')G^{t',t)) 


jk 


(26) 
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From Eq. (26) it is evident that the computational cost scales 
quadratically with the maximum propagation time. The 
quadratic scaling can be reduced to a linear scaling if the col¬ 
lision integral is evaluated in the Markov approximation. The 
Markov approximation consists in neglecting memory effects 
by replacing the pair p(t') and p{t') with the pair p{t) and 
p(t), and in using the equilibrium HE retarded Green’s func¬ 
tion — exp[—i/i^p(f — t')], withthe 

HF single-particle Hamiltonian of the equilibrium system. In 
this case the integral over t' can be done analytically and the 
collision integral becomes a quartic polynomial in p{t). Thus 
the equation of motion for p reduces to a nonlinear differen¬ 
tial equation. In the next sections we benchmark the Markov 
approximation against Configuration Interaction (Cl) and full 
NEGF-i-GKBA simulations. 


IV. HELIUM 

In this section we simulate the TEA spectrum of helium 
atoms recently measured in Ref. 11. Helium is among one of 
the most studied systems in the context of and, 

due to the limited number of electrons, the Cl simulations are 
very accurate. From our point of view helium provides an 
extremely useful platform to benchmark the NEGF-i-GKBA 
methodology. The gas of He atoms is perturbed by a near in¬ 
frared (NIR) transverse pump pulse E(^ t) = 0,0) with 

E{t) = Eq sin^(7rf/Ap) sm{ujpt) for 0 < f < Ap. The ex¬ 
perimental pump intensity is To = 6 x 10^^ W/cm^, which 
corresponds to an electric field Eq = Y^2To/ceo = 6.6 x 10^ 
V/m, the duration of the pump pulse is Ap ~ 15 fs and the 
NIR frequency, ujp ^ 0.57 eV, is slightly detuned from the 
2s — 2p resonance. Thus, the pump alone does not perturb 
the equilibrium state of the He atoms as both the 2 s and 2p 
levels are empty. The situation changes if the probe pulse 
arrives first. In the experiment the probe field is an ultra- 
short pulse e{£t) = (e(t),0,0), with e{t) = eosin^(7r(t — 
t)/ Ap) sm{ujp{t — r)) for r < t < r + Ap. The probe 
pulse has duration Ap ^ 0.5 fs, it is centered at frequency 
ujp = 22 eV and it has an intensity io r\j 10^ W/cm^, which 
corresponds to an electric field eg = 8.6 x 10^ V/m. The den¬ 
sity can been deduced from the pressure P of the gas: 
n(at) = P/{KpT) with ATp the Boltzmann constant. The ex¬ 
perimental estimate of P varies in the range 10 — 240 mbar, 
implying that at room temperature varies in the range 
5.8 X 10^^ — 2.4 X 10^^ cm“^. Finally we observe that the 
experimental thickness (1 pm 1 mm) is much larger than 
the wavelength of the laser pulses (optically thick samples). 
To deal with these thicknesses, we propose the following ap¬ 
proximation dx't) where the ef¬ 

fective thickness ^eff < £ can be used as a fitting parameter. 
As we shall see this approximation well captures the effects 
of screening in the TEA spectrum of a thick helium gas. 

We obtained the one- and two-electron integrals as well 
as the dipole matrix elements from the SMILES package^"^’^^ 
using the VB2 Slater-type orbital (STO) basis, consisting of 
15 basis functions. We performed Cl time-propagations as 
well as NEGF-fGKBA propagations in the HF (S = 0), 2B 
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FIG. 3: (Color online) Density plot of the TEA spectrum (normalized 
to the maximum height) of a helium gas with ^eff = 0.16 mm and 
density = 2.4 x 10^^ cm“^. The pump and probe pulses are 
given in the main text. Results obtained with Cl (top-left), HF (top- 
right), 2B (bottom-left) and the Markovian 2B (bottom-right). 


and Markovian 2B approximation. As a general comment 
we observe that the time-step to achieve convergence in Cl is 
about an order of magnitude smaller than in NEGF-fGKBA, 
thereby Cl requires about ten-times more time steps than 
NEGF-I-GKBA to have the same frequency resolution. 

In Fig. 3 we compare the TEA spectrum in the four dif¬ 
ferent schemes for a gas with ^eff = 0-16 mm and density 
^(at) = 2.4 X 10^^ cm“^. The spectra are obtained by a 
fast Fourier transform of (t) calculated using a time-step 
At = 0.0012 fs and 40.000 time steps in Cl, and At = 0.006 
fs and 7.500 time steps in NEGF-fGKBA. In Cl, see top-left 
panel, the equilibrium peak of the Is — 2p transition occurs at 
frequency 21.1 eV, in agreement with the experiment, whereas 
the Is — 3p transition is not accurate. Nevertheless, as our in¬ 
terest is in the evolution of the Is — 2p transition as the pump- 
probe delay is varied, the VB2 basis is suitable for our pur¬ 
poses. The overall shape of the Cl TEA spectrum is well repro¬ 
duced by all NEGF approximations, see, e.g., the size of the 
splitting and the relaxation toward the equilibrium spectrum. 
The only relevant difference is a uniform frequency shift. We 
also observe that the intensity of the equilibrium peak grows 
monotonically with decreasing r (no coherent oscillations), 
a feature in common with the experiment. In HF, see top- 
right panel, the energy of the Is — 2p transition is overesti¬ 
mated. The inclusion of correlation effects at the 2B level, see 
bottom-left panel, counteracts this overestimation, although 
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FIG. 4: (Color online) Density plot of the TPA spectrum (normalized 
to the maximum height) of a helium gas with ^eff = 1-6 mm and 
density = 2.4 x 10^^ cm“^. The pump and probe pulses are 
given in the main text. Results obtained with Cl (top-left), HF (top- 
right), 2B (bottom-left) and the Markovian 2B (bottom-right). 


the correction is too large. Interestingly the Markovian 2B ap¬ 
proximation, see bottom-right panel, does not change the po¬ 
sition of the HF peak, a result which points to the importance 
of memory effects, see also Section VI. 

The performance of the NEGF-fGKBA approach is sat¬ 
isfactory for larger thicknesses too. In Fig. 4 we compare 
the TPA spectrum for = 1-b mm and density = 
2.4 X 10^^ cm“^ in the four schemes. The fast Fourier trans¬ 
form of has been performed with At = 0.0012 fs 

and 10.000 time steps in Cl, and At = 0.006 fs and 1.500 
time steps in NEGF-i-GKBA. Again all main features of the 
Cl spectrum are well reproduced. The width of the absorption 
peaks is larger as compared to Fig. 3 since the life-time of the 
probe-induced dipole is about an order of magnitude shorter. 

Fet us now come to the physical interpretation of the TPA 
spectrum. For small r the equilibrium peak of the ls—2p tran¬ 
sition undergoes an Autler-Townes (AT) splitting since the 2p 
and 25 levels are mixed by the pump field. The asymmetry 
of the AT intensities is due to the fact that ccp is not exactly 
tuned at the 25 — 2p resonance. For a pump of finite duration 
to generate an AT splitting the probe-induced dipole must de¬ 
cay in a time-window W < Ap. In fact, the pump affects the 
oscillations of dp{t) only in a time-window Ap, hence it can¬ 
not change the position of the peaks of dp{uj) ifW ^ Ap.^^ 
As the time evolution is unitary and the system is finite the 
damping mechanism of the dipole moment is not driven by 


electron-electron scattering. The damping mechanism does 
not have its origin in the radiative decay either (in He relax¬ 
ation through radiative decay is no shorter than hundreds of 
fs). We will explain the origin of the damping mechanism in 
Section IV A. Here we address a different issue, i.e., the bend¬ 
ing of the AT absorption peaks as r decreases and the possible 
formation of a subsplitting structure. 

Consider a three-level He model with basis functions I 5 , 
25, 2p in the oscillating state induced by the probe. If the 
probe is switched off at t = 0 then for t > 0 the probe-induced 
dipole (along the x component) is dp{t) = do sm{uJot), where 
ujQ is the energy of the I 5 — 2p transition and do = dx,is 2 p^^ 
the dipole matrix element. At time — r > 0 we switch on a 
pump field E{t) = sin(co’p(t + r)), of duration 

Ap ^ I/ 7 P, which couples the 2s and 2p levels. We choose 
ccp in resonance with the energy of the 25 — 2p transition and 
work in the rotating wave approximation. Then for times t > 

—T we find^^ dp{t) = (io sin(ccot) cos^-Eotio 

(we assumed for simplicity that the matrix element dx^ 2 s 2 p = 
do). Collecting these results and introducing an exponential 
damping 7 ^ 1/FF (the origin of which is explained in Sec¬ 
tion IV A) we can write 


dpjt) 

do 


= e sin(wot) x 


cos 




t < —T 
t > —T 


(27) 

and dp{t) = 0 for f < 0 (before the probe). This equa¬ 
tion clearly illustrates the behavior previously discussed. For 
t > —T -h I/ 7 P the cosine is essentially constant. Thus, 
if 7 ^ 7 p (hence W Ap) then the pump modifies the 
sin(u;o^) profile only in the time window (— r, — r + 1/7p), 
too short to change the position of the peaks in Accq of the 
Fourier transform dp{uj). Fet us now consider the opposite 
limit 7 ^ 7 p (hence W <C Ap). For all times t smaller 
than 1/7 (after this time the dipole is exponentially small) we 
can approximate the cosine with cos{Eodo{t + r)). In this ap¬ 
proximation the Fourier transform dp{uj) has a simple analytic 
form and for, e.g., uj ujo E Eodo, we find 


7 e 7cos(cjr) + (cj - Eodo) sin(cjr) 

-ImK(„)] == - - + - 

(28) 

with uj = cc—cco- The denominator of Eq. (28) has a minimum 
in 00 = Eodo which is independent of r. However, the maxi¬ 
mum of —lm[dp{uj)] does not coincide with the minimum of 
the denominator as r decreases from zero. For small negative 
r the maximum occurs at frequencies (j Eodo{lErj/2) < 
Eodo. This explains the bending of the AT absorption peaks (a 
similar analysis can be done for frequencies cc ujo — Eodo). 

Noteworthy, the full Fourier transform of the dipole in 
Eq. (27) yields spectra that resemble very closely those in 
Figs. 3 and 4. In Fig. 5 we display the density plot of 
—ujlm[dp{uj)] (X &{uj) for cco = 21.3 eV, Eodo = 0.5 eV, 
7 p = 0.41 eV (corresponding to a duration Ap ^ 10 fs) 
and for 7 = 0.21 eV (left panel, compare with Fig. 3) and 
7 = 1.37 eV (right panel, compare with Fig. 4). Equation (27) 
does therefore provide a convenient analytic parametrization 
of experimental TPA spectra. For longer pumps Eq. (27) pre- 
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FIG. 5: (Color online) Density plot of —ulm[dp(uj)]/dQ, with dp the 
Fourier transform of the function dp in Eq. (27), for ujq = 21.3 eV, 
Eodo = 0.5 eV, 7 p = 0.41 eV and for 7 = 0.21 eV (left panel) and 
7 = 1.37 eV (right panel). 


diets the formation of a subsplitting structure as well. In Fig. 6 
we show the TPA spectrum as a function of the dipole life¬ 
time 1/7 at delay r = 0 for ccq = 21.3 eV, Eodo = 0.2 
eV, 7 p = 0.125 eV (corresponding to a duration Ap ^ 33 
fs). In addition to the AT peaks at ccq i Eodo two extra peaks 
emerge, in agreement with recent experiments on a thick he¬ 
lium gas perturbed by a NIR pump of duration 33 fs.^^ We 
point out that the number of extra peaks increases with in¬ 
creasing the AT splitting (hence with increasing the intensity 
of the NIR pulse), a prediction which could be easily checked 
experimentally. 


A. Damping in closed systems with unitary evolution 

The damping of the probe-induced dipole in the helium gas 
has been systematically investigated both numerically and ex¬ 
perimentally in Ref. 66 . Here we provide a transparent ex¬ 
planation based on the analytic solution of the Schrodinger 
equation ( 8 ) for a simple two-level model with two-particle 
states |ls^) and \ls2px). We take the equilibrium Hamilto¬ 
nian H diagonal on this basis and let ei and 62 be the corre¬ 
sponding eigenenergies. We denote by do = {ls‘^\dx\ls2px) 
the x-component of dipole matrix element and write the 
dressed probe field along x in accordance with Eq. (11), i.e., 
Ep{t) = e{t) -h ad{t), where a = > 0 

d{t) = {t)\dx\'^{t)). Expanding the time-dependent two- 

particle state as = ci(t)|ls^) + C 2 {t)\ls 2 px) we find 

^ do£p{t)c 2 {t), (29a) 

ic 2 {t) = e 2 C 2 {t) ^ do£p{t)ci{t). (29b) 

Eor any real £p the time evolution of the coefficients ci and 
C 2 is unitary. The expression of the time-dependent dipole 
in terms of ci and C 2 reads d{t) = 2doRe[cl{t)c2{t)]. We 
can get a differential equation for d if we introduce two more 
real functions g{t) = Im[cJ(t)c 2 (t)] and f{t) = |ci(t)p — 
|c 2 (t)p. It is a matter of simple algebra to show that after 
the probe (hence e{t) = 0 ) these three functions satisfy the 



FIG. 6 : (Color online) 3D plot of —ulm[dp{uo)] /do (normalized to 
the maximum height), with dp the Fourier transform of the function 
dp in Eq. (27), as a function of the dipole life time I /7 at delay r = 0 
for cjo = 21.3 eV, = 0.2 eV, 7 p = 0.125 eV. Arrows indicate 
the extra peaks discussed in the main text. 


system 


f{t) = -2doag{t)d{t), (30a) 

^ (30b) 

zdo 

d{t) = 2doujogit), (30c) 

where we introduced ccq = — 62 . Taking the time derivative 

of Eq. (30c) and using Eq. (30b) we find 

d{t) + uJod{t) — 2ad^uJof{t)d{t) = 0. (31) 


With the help of Eq. (30c) we rewrite Eq. (30a) as f{t) = 
— ^d?{t). Therefore /(t) is a monotonically decreasing 
function of time and, by definition, it is bounded between —1 
and 1 . This implies that limt^oo f{t) = f00 C (-1,1). If /oo 
were positive then the long-time solution of Eq. (31) would be 
an oscillatory function with an exponentially increasing am¬ 
plitude, in contradiction with the fact that d{t) G {do, —do). 
We conclude that the limiting value foo G (—1,0) indepen¬ 
dently of the initial condition. Consequently, for large t the 
function d{t) oscillates at frequency ccq with an amplitude de¬ 
caying as where 7 = 2 Q((io<^o|/oo|- Otir analysis does 
explain the physical origin of the damping as well as the de¬ 
pendence of 7 on the thickness and density of the gas. In fact, 
7 oc a (X therefore the larger the thickness and/or 

the density is and the faster the amplitude of the dipole oscil¬ 
lations decays. 


V. PUMP-INDUCED IONIZATION 

In this Section we extend the NEGE+GKBA formalism to 
deal with ionization processes induced by the pump. Eor this 
purpose it is convenient to work with the HE orbitals. Let 
be the equilibrium density matrix and 

~ Prim (^^) 


mn 
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the equilibrium HF Hamiltonian in the original basis. The HF 
orbitals a^ipi{Y) diagonalize both and 

and are orthonormal. To distinguish the HF basis from the 
original basis we use greek letters to label the former. We 
have 


^HF,/xzy ~ (33) 

and = 5 where = 1 if < ep < 0 and = 0 
otherwise, ep being the Fermi energy. 

In finite systems like atoms and molecules the HF orbitals 
with > 0 are states in the continuum. These are the states 
that get occupied by the photoelectron in a ionization process. 
We assume that the Coulomb interaction between photoelec¬ 
trons and bound electrons is negligible and set to zero the two- 
electron integrals with at least one of the four indices 

in the continuum (this amounts to neglect Auger transitions). 
Then the self-energy has nonvanishing matrix elements only 
between bound states. In Appendix A we prove that the equa¬ 
tion of motion for the density matrix with both indices in the 
bound sector reads 

+ [hHFit), p{t)] = i [I{t) + Iion{t)] - H.C.. (34) 

In Eq. (34) every matrix has indices running over the bound 
states. The integral /ion accounts for the pump-induced ion¬ 
ization (since Tr[/ion + H.c.] / 0 the number of bound elec¬ 
trons is not conserved) and it is calculated like in Eq. (23) 
except that the correlation self-energy is replaced by the ion¬ 
ization self-energy Sion- The latter has a vanishing lesser part 
and a greater part given by (see Appendix A) 

(35) 

Here £i is the i-th component of the electric field S = 
{£x,£y, £z) and the tensor 

- i') = “* S ) dj^au, (36) 

aGc 

where di is the i-th component of the vector of matrices d = 
(dx^dy^dz) and the sum over a runs in the continuum. In 
Fourier space 


O',] 


,(w) = -2TTi di^^a S{ui - €a) dj^, 

1 


aEc 


2i E 


aEc 


CC - Cq, -h if] 


dj,au^ (37) 


where is a positive constant of the order of the level spacing 
of the continuum states. Typically the ionization is caused by 
the action of a pump pulse with a Fourier transform peaked 
around some frequency ujp (larger than the ionization energy 
of the system). Therefore, E?^^ is dominated by those terms 
in cr(t — t') that oscillate at frequency Cq, ujp. By virtue of 
this observation we implement a time-local approximation 

~ ( 38 ) 


which implies — t') = cr'^^^{ujp)5{t — t'). Substituting 

this result into Eq. (35) yields 

= (39) 


where 


( 40 ) 

is a self-adjoint positive-definite matrix for all times t. 

The time-local approximation allows us to simplify the in¬ 
tegral /ion(^) appearing in Eq. (34). Taking into account that 
E^^n = 0 we have 

Iion{t)= f dt'^>,{t,t')G<{t',t) 

J — oo 

= r{t)p{t). (41) 


Inserting this result into Eq. (34) we finally obtain 


+ [^HF(i), p{t)] - i {r(t), = il{t) - H.c. 


(42) 

where the curly brackets signify an anticommutator. Equa¬ 
tion (42) constitutes the generalization of the NEGF-fGKBA 
formalism to open systems. 


VI. KRYPTON 

We apply the formalism of the previous Section to address a 
retardation effect observed by Goulielmakis et al.^^ in the TPA 
spectrum of a krypton gas. In the experiment a strong pump is 
shone on the gas, electrons from the 4p shell are expelled and 
an admixture of Kr atoms and ions, with n = 1, 2, 3 ..., 

is formed. The admixture is subsequently probed with an ul¬ 
trafast pulse, thus inducing transitions from the 3d to the 4p 
shell. The main focus in Ref. 34 was on the coherent oscilla¬ 
tions of the peak intensities of the Kr^+ ion as a function 

of the pump-probe delay r. However, the experimental TPA 
spectrum reveals another interesting feature as a function of 
r. The absorption peaks of Kr^+ develops after the absorp¬ 
tion peaks of Kr^+, implying that it is faster to expel one elec¬ 
tron than two electrons. To reproduce this retardation effect 
theoretically a formalism for the TPA spectrum of an evolv¬ 
ing admixture is needed. The NEGF-fGKBA approach is in 
principle suitable for this purpose. As we shall see, important 
qualitative aspects of the admixture are intimately related to 
the diagrammatic structure of the self-energy. 

The Kr gas is ionized by a few-cycle NIR pump E(/t) = 
(/^(t),0,0) with E{t) = /^o sin^(7rt/Ap) sin(cc;pt) for 0 < 
t < Ap. The experimental pump intensity is To = 7 x 10^^ 
W/cm^, corresponding to an electric field Et^ = = 

7.2 X 10^^ V/m, the duration of the pump pulse is Ap ^ 
7.6 fs and the NIR frequency is ccp ~ 1.65 eV. After a 
time r the Kr admixture is probed with an extreme ultra¬ 
violet attosecond pulse e(/t) = (e(t),0,0), with e(t) = 
eosin^(7r(t — f)/Ap) sin(co’p(t — r)) for r < t < f -b A^. 



9 



FIG. 7: (Color online) Transient ionization yield per spin (solid line) 
in HF and 2B (indistinguishable) and amplitude of the pump pulse 
(dashed line) in arbitrary units. 


Here r = f — (Ap — Ap)/2 is the time-distance between the 
maxima of the pump and probe pulses. The probe pulse has 
duration Ap ^ 150 as, it is centered at frequency ujp = 80 eV 
and it has an intensity io ^ 10^^ W/cm^, which corresponds 
to an electric field eo = 8.6 x 10^ V/m. We discard the dress¬ 
ing of the probe field and solve the equation of motion for p 
with S = e. 

The one- and two-electron integrals as well as the dipole 
matrix elements have been calculated with the SMILES pack- 
age^4,65 using 66 STO basis functions taken from Ref. 67. 
As we are not interested in the coherent oscillations of the 
peak intensities we do not include the spin-orbit coupling re¬ 
sponsible for the splitting of the 4p and 3(i orbitals. Thus 
we should expect one main absorption peak per ion, corre¬ 
sponding to transitions from the Ap to the 3(i shell. We find 
eighteen HF states with energy below zero. The remaining 
HF states are used to construct the ionization self-energy ac¬ 
cording to Eq. (37). The simulations show that electrons are 
essentially removed from the Ap shell in agreement with the 
analysis of Ref. 34. In Fig. 7 we display the transient ioniza¬ 
tion yield, i.e., the expelled charge per spin, during the action 
of the pump. The charge is expelled in pockets at a rate of 
twice the frequency of the laser pulse, in agreement with the 
Cl calculations. ^^’^4 Interestingly, the HF and 2B yields are 
indistinguishable. The situation is drastically different for the 
TPA spectrum, with the 2B approximation performing much 
better than the HF one (see below). Whether the absorption 
onset of Kr^+ is earlier than the absorption onset of Kr^+ is 
the central issue addressed below. 

In the upper panel of Fig. 8 we report the TPA spectrum of 
the Kr admixture in the HF approximation. We have propa¬ 
gated the system for ^ 50 fs after the action of the probe with 
a time step At = 0.0025 fs (corresponding to 20.000 time 
steps) and then broadened the Fourier transform of dp{t) by 
0.8 eV to account for the experimental resolution. The result 
is extremely disappointing. The HF TPA spectrum constitutes 
of one doublet (merged in Fig. 8 in one single line-shape due 
to the broadening) with a simultaneous raise of both peaks. 
The peaks correspond to transitions from the 3d shell to ei- 



CO (eV) 92 



FIG. 8: (Color online) TPA spectrum (normalized to the maximum 
height) of a krypton gas in the HF (top panel) and 2B (bottom panel) 
approximation. White (dashed) lines in the bottom panel are guide 
for the eyes to better visualize the retardation effect. The pump and 
probe pulses are given in the main text. 


ther the Ap^ orbital or the Apy^^ orbitals. In fact, these tran¬ 
sitions are nondegenerate in the HF approximation. As the 
pump is polarized along x the Ap^ orbital looses more charge 
than the Apy^z orbitals, thereby breaking the degeneracy (al¬ 
beit only slightly). Even more noteworthy, however, is the 
absence of spectral structures due absorption of multiply ion¬ 
ized Kr atoms. The numerical simulation has been repeated 
with pump fields of different frequencies and intensities but 
no sign of other absorption peaks has been observed. The first 
important conclusion of this preliminary study is that the ap¬ 
pearance of absorption peaks in multiply ionized Kr atoms is 
a correlation ejfect. 

We have then included correlation effects at the level of 
the Markovian 2B approximation but the outcome has not 
changed (not shown). The main difference between the HF 
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and the Markovian 2B spectra is an overall frequency shift. 
The so far accumulated numerical evidence leads us to con¬ 
jecture that any time-local approximation (no memory) to the 
self-energy is doomed to fail. We mention that the equation of 
motion for p with a time-local S has the same mathematical 
structure of the Time-Dependent Density Functional Theory 
(TDDFT) equations at the level of the Adiabatic Local Den¬ 
sity Approximation (ALDA). Hence, TDDFT spectra at the 
ALDA level would also fail in capturing the absorption peaks 
of multiply ionized atoms. The second important conclusion 
is that static correlation effects are not enough. 

Dynamical correlation effects are contained in the full (non¬ 
local in time) 2B self-energy. We have solved the equations of 
motion for p with the collision integral of Eq. (26). The TPA 
spectrum is shown in the lower panel of Fig. 8. We clearly dis¬ 
tinguish two structures corresponding to the TPA spectrum of 
Kr^+ and Kr^+. In fact, the energy gap between the structures 
is consistent with the ~ 3 eV experimental gap of these two 
ions.^^ Remarkably, the high-energy structure develops ^ 3 
fs after the low-energy one. This is the aforementioned retar¬ 
dation effect which we have just proved to be within reach of 
the NEGF-fGKBA approach. Furthermore, the obtained delay 
is commensurate with the 5 fs delay observed in experiments. 
A non-local in time self-energy is crucial for the appropriate 
description of the Kr (and probably of any other) evolving ad¬ 
mixture. 

Although the 2B approximation represents a noticeable im¬ 
provement over time-local approximations there still remains 
one issue to address. The experimental TPA spectrum con¬ 
tains small absorption peaks attributable to transitions in the 
Kr^+ ion. We have not been able to see these structure within 
the 2B approximation. Although we are not aware of any for¬ 
mal result relating the possibility of describing multiply ion¬ 
ized atoms to the diagrammatic structure of the self-energy we 
observe that the kernel 6'E/6G contains at most one particle- 
hole excitation in HE and two particle-hole excitations in 2B. 
It is therefore tempting to argue that in order to observe the 
absorption peaks of Kr’^+ the kernel 6T^/6G should contain 
at least n particle-hole excitations,^^’^^ which implies that S 
should contain diagrams of order at least n in the interaction. 
Finding a general solution to this problem would certainly 
be valuable and contribute to advance the understanding of 
many-body diagrammatic theories. 


VII. CONCLUSIONS 

We have introduced a NEGF-fGKBA approach to transient 
photoabsorption experiment suitable for pump fields of arbi¬ 
trary strength, frequency and duration and for any delay be¬ 
tween pump ad probe pulses (hence for delays in the overlap¬ 
ping regime too). The size of the arrays in NEGF calculations 
scales quadratically with the number of basis functions. The 
Coulomb interaction between electrons is included diagram- 
matically through the correlation self-energy and the possi¬ 
bility of ionization can be described through the ionization 
self-energy. 

The approach has been benchmarked against the TPA spec¬ 


trum of He reported in Ref. 11. Helium is a weakly correlated 
system and all self-energy approximations have been shown to 
agree with the Cl results. We have provided a simple yet rig¬ 
orous explanation of the bending of the AT absorption peaks 
and derived a useful formula for fitting the experimental TPA 
spectra. We have also addressed the exponential damping of 
the probe-induced dipole and related it to the thickness and 
density of the gas. 

A more severe test for the NEGF-i-GKBA approach is the 
TPA spectrum of Kr reported in Ref. 34. We have shown that 
for a proper description of the evolving admixture of the Kr 
ions the self-energy should have memory. This is not the case 
for the HE and Markovian 2B self-energies which yield the 
TPA spectrum of a pure Kr^+ ion. We argue that the situation 
does not change in TDDFT with ALDA exchange-correlation 
potentials. On the contrary, the full 2B self-energy leads to a 
second structure in the TPA spectrum that is assigned to Kr^+ 
and that develops about 2-3 fs after the first, in fair agreement 
with the experiment. More theoretical and numerical work 
is needed to understand the relation between self-energy dia¬ 
grams and the emergence of absorption peaks due to multiply 
ionized atoms. 
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Appendix A: The embedded GKBA equation for p 

The lesser and greater Green’s functions follow from the 
Keldysh Green’s function G{z^z') with arguments z and 
z' on the Keldysh contour. In particular G^ (G^) is the 
Keldysh G with the first (second) contour argument on the 
forward branch and the second (first) contour argument on the 
backward branch. The Keldysh G satisfies the equations of 
motion^ ^ (in matrix form) 

[i-^ - hHF{z)]G{z,z') = 6{z,z') 

+ JdzT,{z,z)G{z,z'), (Al) 

G{z,z')[-i-^ - hiiF{z')] = S{z,z') 

+ JdzG{z,z}T,{z,z'), (A2) 

where the integral is over the Keldysh contour. Choosing z 
on the backward branch and z' on the forward branch, and 
applying the Langreth rules we obtain the equations of motion 
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forG< 

+ jdt'E<{t,t)G^{t,t'). (A3) 


[- - ^•HF(i')] = JdiG^{t,t)T,<{t,t') 

+ J diG<{t,t)T,^{t,t'). (A4) 

The upper indices “R” and “A” signify retarded and advanced 
functions respectively. These are defined according to 

= ±e{±tTt')[F>{t,t') - F<{t,t')]. (A5) 

Subtracting Eq. (A4) from Eq. (A3) and setting t' = t we 
find the equation of motion Eq. (20) for the density matrix 
p{t) = -iG<{t,t). 

Let us work in the HE basis of the equilibrium system 
and write + 5p^y{t). In the same basis the 

HE Hamiltonian in Eq. (21) reads 


We can make use of the block structure of h^F and E to 
simplify the equations of motion for the Keldysh G. In the 
bound-bound sector Eq. (Al) reads 


- h’^!p{z)]G’>\z,z') - [£{z) ■ d^^]G^\z,z') 

= S{z,z') + JdzE'^\z,z)G^\z,z'), (All) 

whereas in the continuum-bound sector the same equation 
reads 


- h'^p{z)]G^\z,z') - [£(z) ■ d-^]G^\z,z') = 0. 

(A12) 

We define the continuum (noninteracting) Green’s function 
as the solution of 


[i- - z') = S{z, z'), (A13) 

and rewrite Eq. (A 12) in integral form 

G'^\z, z') = j dz g'^^z, z) [£{z) ■ d"'’] G'^\z, z'). (A14) 


^HF,/xzy(^) = ' d^zy (A6) 

with 

flFLF 

aP 

“ 1 “ ^ ^ PPajt^' (^^) 

a/3 


Inserting Eq. (A 14) into Eq. (All) we find 

[i±-h1^^{z)]G^\z,z')=5{z,z') 

+ Jdz [i:’>\z, z) + SfUz, I)] , z'), (A15) 

with the ionization self-energy defined according to 


The HE states can be grouped according to their energies: if 
< 0 then is a bound state, otherwise is a continuum 
state. We assume that the two-electron integrals v^ocpu with at 
least one index in the continuum are negligible and set them 
to zero. Consequently, w^apu with at least one index in the 
continuum vanishes too. Let us represent a matrix M. in the 
HE basis as 


M = 




(A8) 


where in both indices run over the bound states, in 
the first index run over the bound states and the second index 
over the continuum states, and so on. Then, the HE Hamilto¬ 
nian has following block structure 


Hflf 


£ • d'>'= \ 

£ ■ d^’’ (ifp ) ’ 


(A9) 


where we took into account that h^p = 0, see Eq. (A7). Sim- 
ilarly, we infer that the block structure of the correlation self¬ 
energy is 


^’^Uz,z') = [£{z)-d^’^Y%z,z')[£{z')-d% (A16) 

Thus, the continuum states can be downfolded in an exact way 
into an effective equation for . A similar equation can be 
derived starting from Eq. (A2) and reads 

G'’\z,z')[- - h%,{z')] = 5{z,z') 

IdzG^\z, z) [y!>\z, z') + E\l{z, ^0] • (A17) 

Below we use Eqs. (A15, A17) to generate an equation for 
the density matrix in the bound-bound sector. To lighten the 
notation we omit the upper indices “66”, so a matrix with no 
upper indices is a matrix in the bound-bound sector. 

Comparing Eqs. (A15, A17) with Eqs. (Al, A2) we deduce 
that the equations of motion for G^ are the same as Eqs. (A3, 
A4) except that the correlation self-energy is replaced by E + 
Sion- Therefore, the equation of motion for p is the same as 
Eq. (20) except that the collision integral is calculated with 
Erom Eq. (A13) we have 


E = 


aGc 


s'*'’ 0 

0 0 


(AlO) 


(A18) 
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where W^{t) = T[e“*>^o the evolution operator 

in the continuum sector whereas and = 1 — nc,. 

For a in the continuum we have Cq, > 0 > ep and hence = 
0 (which implies = 0) and = 1. The evolution 

operator takes a very simple form if we ignore the effect of 
the pump between continuum states, i.e., if we approximate 
0. In this case see Eqs. (A6, A7), and 

hence 

= (A19) 


From Eq. (A 16) it follows that the greater ionization self¬ 
energy is 




a. 

which agrees with Eqs. (35, 36). 
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